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INTRODUCTION 


Most  of  our  results  have  been  documented  in  detail  in  the  indi¬ 
vidual  progress  and  annual  reports.  Our  research  and  activities  have 
included  direct  studies  in  support  of  gravity  gradiometer  survey  missions 
as  well  as  participation  in  related  activities  in  wh-'ch  gravity  mapping, 
instrument  development  work,  reviews  of  gravity  gradient  programs  and 
technical  conferences  involving  the  development  of  gravity  gathering 
equipments,  and  reviews  of  states  of  the  art.  Indirectly  related  instru¬ 
mental  work  in  gravity  gradiometry  has  included  cryogenic  gravity  gradi¬ 
ometer  development  supported  from  AFOSR  and  jointly  funded  with  our 
Physics  Department.  This  class  of  instrument  and  its  related  technology 
is  a  candidate  both  for  future  gradiometers  for  survey  work  in  aircraft 
and/or  in  orbit  as  well  as  for  future,  more  accurate  navigation  require¬ 
ments.  We  have  benetited  from  studies  we  have  performed  for  Goddard  Spac 


Flight  Center  and  the  Johns  Hopkins  Applied  Physics  Laboratory  in 
connection  with  GRAVSAT  and  gradiometers  for  orbital  measurements. 

In  addition  to  our  direct  support  for  their  work,  we  have  participated 
on  planning  groups  and  chaired  a  workshop  at  Goddard  reviewing  the 
technologies,  both  room  temperature  and  low  temperature,  for  orbital 
gravity  gradiometry.  We  have  participated  in  a  review  of  the  National 
Academy  of  Engineering  of  Precise  Navigation  for  the  Navy.  Included 
in  this  advisory  function  was  a  study  of  the  use  of  low-temperature 
technology  for  future  Navy  needs.  Gravity  gradiometers  played  an 
important  role  in  all  future  navigation  system  considerations.  In 
addition,  we  participated  in  some  internal  DOD  reviews  of  technology 
strategies  to  ensure  the  best  possible  results  for  investment  in 
development  in  gradiometry. 


There  are  a  variety  of  issues  associated  with  data  gathering  and 
fitting  gravity  fields  which  have  been  studied.  The  aliasing  due  to 
spacing  between  ground  tracks  is  a  subject  appropriate  both  for  satel¬ 
lite  data  gathering  as  well  as  low  altitude  mapping  missions.  The 
local  weighting  of  GRAVSAT  doppler  data  in  estimating  surface  gravity 
anomalies  and  undulations  were  reported  on.  The  methods  employed  have 
application  to  other  gravity  gathering  missions.  An  important  step  in 
changing  the  methods  of  estimating  accuracy  of  gravity  model  fitting 
was  made  by  inverting  the  classical  problem.  Typically,  earlier  studies 
were  done  by  doing  sensitivity  analyses.  Varying  a  parameter  such  as  a 
spherical  harmonic  coefficient  of  the  gravity  field  and  looking  at  the 
effect  on  a  measurement  determined  that  there  would  be  adequate  sensi¬ 
tivity  if  a  response  was  larger  than  the  noise.  However,  many  parameters 


exist  in  the  gravity  field  and  this  does  not  indicate  whether  they  can 
be  separated  or  not.  A  study  on  estimating  accuracy  of  gravity  gradient 
survey  systems  developed  a  theory  of  two-dimensional  Fourier  representa¬ 
tion  which  changed  the  problem  to  the  more  direct  estimation  theory. 

Here  one  can  answer  questions  of  observability  and  determine  directly 
the  accuracy  with  which  one  could  expect  each  of  the  parameters  to 
be  measured. 

A  variety  of  parameterizations  are  possible  with  gravity  fields. 
Classically,  spherical  harmonic  expansions  have  been  used  because  the 
expansions  involve  orthogonal  functions.  The  advantage  here  is  that 
one  can  truncate  such  theories  and  know  that  the  remaining  functions 
represent  the  best  fit.  By  contrast,  the  spherical  harmonic  expansion 
involves  the  gravity  around  the  entire  sphere.  For  much  of  the  current 
interest  in  gravity,  small  patches  must  be  modeled.  For  high  resolution 
the  number  of  parameters  needed  becomes  excessive  when  a  world  wide 
parameterization  is  used.  One  solution  to  this  is  to  use  mass  con¬ 
centrations  and  the  number,  size,  and  location  of  masses  then  become 
the  parameters  with  which  one  works.  Locally  the  size  and  variety  is 
intense  whereas  at  greater  distances  they  can  be  sparse  and  still 
provide  the  resolution  desired.  An  earlier  report,  "MASCON 
APPROXIMATIONS,"  has  been  revised  during  the  past  year  (August  1985) 
and  is  included  as  an  enclosure  with  this  report. 

Designing  a  gravity  survey  can  have  a  significant  influence  both 
on  the  cost  and  the  accuracy  of  the  resulting  model.  Many  investigators 
have  made  contributions  in  terms  of  spacing  and  the  proper  utilization 


of  data  gathered  in  order  to  make  the  numerical  data  reduction  tech¬ 


niques  efficient.  During  the  closing  portions  of  this  contract,  we 
have  addressed  the  question  of  the  pattern  that  would  provide  the  best 
opportunity  to  separate  instrument  noise  from  uncertainties  in  the 
gravity  field.  If  one  takes  an  orthogonal  set  of  tracks  and  keeps 
that  pattern  constant,  the  question  arises  whether  the  order  in  which 
the  tracks  are  traced  would  make  a  difference  in  the  separability  of 
instrument  noise  from  the  gravity  uncertainties.  Track  crossings 
are  when  the  measurements  from  the  instrument  should  be  the  same  along 
each  path.  The  conjecture  was  that  it  would  since  the  crossings  occur 
with  a  different  distribution  in  time.  Thus,  a  pattern  in  which  no 
crossings  occurred  until  late  in  the  survey  would  leave  no  check  points 
during  the  first  half  of  the  survey,  for  example.  By  contrast,  skipping 
paths  would  allow  one  to  get  the  first  crossing  much  earlier  and  dis¬ 
tribute  the  crossings  more  uniformly  throughout  the  survey.  A  theory 
was  developed  to  evaluate  the  relative  expected  values  of  error  and 
some  numerical  examples  were  calculated  to  indicate  the  promise  or 
lack  thereof  of  success.  Unfortunately,  the  results  do  not  indicate 
the  improvement  would  be  significant.  This  theoretical  development 
and  the  numerical  examples  have  been  prepared  in  a  paper  which  will 
be  presented  at  the  Air  Force  Academy  on  February  11,  1986.  A  copy 


of  the  paper  is  enclosed  as  a  portion  of  this  report. 
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MASCON  APPROXIMATIONS 


John  V.  Breakwell 
Department  of 
Aeronautics/  Astronautics 
Stanford  University 
Stanford,  CA  94305 


by 

Weilian  Yang 
Chinese  Academy  of 
Space  Technology 
Beijing 

Peoples  Republic  of  China 


I.  INTRODUCTION 

The  error  in  representing  the  earth’s  gravity  field  by  a  large  but  finite 
number  of  mascons  is  clearly  concentrated  in  the  short  wavelength  part  of  the 
gravity  spectrum.  It  is  thus  appropriate  to  use  “flat-earth”  approximations,  as  in 
[1],  to  estimate  this  error. 

The  most  appropriate  thin  layers  in  which  to  locate  the  mascons  are  the 
idealized  oblate  spheroidal  surface  and  the  Haddon-Bullen  layer  15  km  below  the 
surface,  at  the  bottom  of  the  earth’s  crust.  Assuming  that,  as  in  [2],  the  earth 
density  follows  Jeffreys’  model  of  regional  isostatic  compensation,  the  surface 
density  fluctuations,  confined  to  these  two  layers,  can  be  estimated  as  in  [l]  from 
data  of  various  sorts  obtained  from  low  satellites  in  polar  orbit.  .Assuming  a  grid 
of  mascons  in  each  layer  with  separation  A  both  .VS  and  EW,  the  mass  of  a 
typical  mascon  may  be  chosen  as  the  estimate  of  the  integral  of  the  surface 
density  over  a  square  of  side  A  centered  at  the  mascon  grid  point. 


It  is  shown  in  this  paper  how  to  calculate  the  two-dimensional  spectrum  of 
the  additional  error  in  such  quantities  as  the  vertical  and  horizontal  gravity 
components  at  satellite  altitude  z  due  to  the  finite  separation  A  of  the 


mascons. 
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Numerical  results  are  presented  for  the  case  where  the  density  fluctuations 
are  estimated  from  gradiometer  measurements  taken  in  a  low  satellite,  evaluated 
for  various  ratios  A  j z. 


2.  THE  2-LAYER  MASCON  MODEL 

According  to  Jeffrey’s  dynamical  model,  the  disturbing  potential  can  be 
expressed  as 


u(r,y,r)  =  7  /  / 


pt(u,v)dudv 


+Kf  / 


[(*-») 2  +  (SM)2  +  ;2]  1/2 
Pn(  U,l')dudv 


(2.1) 


[(z-u)2  -I-  (y-v)  2  +  (r+tf)  2]  1/2 

where  px[u,v)  is  the  surface-density  fluctuation  of  the  surface  of  the  earth  and 
P2(u,v)  is  the  surface-density  fluctuation  of  a  compensating  layer  at  depth 
H=  15  km  [2].  7  is  the  universal  gravitation  constant.  The  relation  between 
p{(u, r)  and  p«{u,v)  is 


d  2 

d-  )  2 

du  2 

+  a,2  1  +IJ 

pz{u,r)  =  -/^(u.r) 


o  o\ 


where  d  =  48.1  km  in  accordance  with  the  Haddon-Bullen  earth  model,  (2.1)  can 
be  rewritten  as 


OO  00 


u(i,y,z)  —  7  £  E  /  J 


pl(u,v)dudv 


m— oo  D.m  [(*-«)  2  +  (V-V)  2  +  Z2}  1/2 

p2(  u,v)dudv 


(2.3) 


OO  OO 


+  i  E  Elf 


n— -oo  m»-oo 


/>..  ((*-«) 2  +  {y-v) 2  +  (*+#) 2] 1/2 


where  the  integral  areas  Z)nm  are  defined  as 

m  =  {(x,y)|(n  -  l/2)Ax  <  x  <  (n  +  1/2 )Ar,  (m  -  1/2)A y  <  y  <  (m  +  l/2)Ay} 

(2.4) 

If  Ax, Ay  are  small  enough,  we  can  approximate 

pl(u,v)dudv 


bv 


d!  [(*-«) 2  +  (y-t>)2  +  z2]  1/2 
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and 
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D.m  ((x-u)2  +  (y-v)2  +  (;+//) 2]  1/2 


by 


A/ (2) 
iVlnm 


{ [x  -  nAx]  2  +  [y  -  mAy]  2  +  (z  +  H)  2}  '/2 
and  obtain  the  potential  a^(x,y,x)  of  the  2-layer  mascon  model: 


(2.5) 
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where 


& 
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where  //,  =  0,  //2  =  H.  The  simplest  choice  of  A/j^  is  AxAjrtpj?),,,*,  where 
(p\%nm  ’s  the  estimate,  from  satellite  tracking  or  otherwise,  of  the  average  over 
Dnm  of  surface  density  pt.  We  shall  be  mainly  concerned  in  this  paper  with  this 
definition  of  A/J^,  hut  in  the  final  section  we  shall  look  at  more  general 
definitions  of  average  density. 

We  may  express  u(x,y,.:)  and  u*{x,y,;)  as  two-dimensional  convolutions: 


u{z,y,z) 


=  V 


<1  yr'+  J”  +(.’+  //,)  ' 


ul(x,y,z) 


2 

=  V 


«-i  v/j‘ +  y  +  u  +  Hx) 


=■  ®  (/T,  (z,y)£>(x.j/)|  ,  (2.8) 


where  (x.y)  denotes  the  estimate  of  the  average  density  in  layer  i  over  a 
rectangle  with  center  (x.y)  and  dimension  Ax  Ay,  and  D[x,y)  denotes  the 
periodic  function: 


D(x.y)  =  AxA y  V  V]  S(x  -  nAx)<5(y-  mAy) 

-oo  m— -oo 


with  period  Ax  in  x  and  period  A y  in  y. 

This  periodic  function  is  easily  found  to  have  the  Fourier  series 
representation: 


oo  oo  2ff;|  —  +  - 

D{x, y)  ~  V  V  e 


(2.10) 
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3.  ERROR  ANALYSIS 


Assume  that  pt{x,y)  (t  =  1,2)  are  stationary  random  processes.  Using 
two-dimensional  Fourier  transformation  as  in  reference  [1],  (2.7)  becomes 

=  —  «'2W[PiP)  +  «-%2P)]  ,  (3.1) 

(jJ 


where  3  =  ,  u  = 


+  w  2  and  pp)  denotes 


j7  pl(i,y)e";lu'r,  +  "’,')dxdy 


The  dynamic  compensating  relation  (2.2)  becomes: 


po(Zj)  — 


~  PiP) 

1  +  <f4w4 


Combining  (3.1),  with  z  =  0,  and  (3.2),  we  obtain 


/>P)  =  Cp)U0p)  , 


where  t '0(3J)  denotes  (p, 0),  and 


C,P)  = 


2?T7  1 


C,(3)  = 


The  mascon  model  potential  becomes: 


27n[e  HjJ  -  (l-Hf4^4)] 


^(3,:)  =  e-*'[plHZ3)  +  e-u'p2*(ZJ)\  , 


$ 
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where  />,*p)  denotes  the  Fourier  tranform  of  [/T,-  (z,y)D(z,y)],  expressible 
because  of  (2.10)  as 


where 


00  OO 


pfa)  =  e  e  r,-(3-aj  , 


n— -oo  m^-oo 


uJ  = 

wnm 


The  aliasing  of  p  ,•  in  (3.6)  is  due  to  the  finite  mascon  grid  size. 

In  general,  some  gravity-related  quantity  q(W),  in  which  we  are  interested, 
can  be  expressed  as 

»lP)  =  /,PV,P)  +  /2P)/>2P)  (3.8) 

or,  using  matrix  notation, 

n(2)  =  /rP)pP)  ,  (3.9) 

where 

/lP)  /»i(0) 

TP)  =  ,  p(  3)  =  (3.10) 

/2P)  P2P) 

And  the  approximation  77  *p),  of  r/(3)  in  the  mascon  model  is  given  by 

v*P)  =  /rP)  E  E  P  (Hm)  •  (3.11) 

n— -oo  m— -oo 

The  next  step  is  to  analyze  the  error  At]  caused  by  using  the  mascon 
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model,  where 


A//(ZD)  =  r}  t(Z3)  -  »/(aJ) 


(3.12) 


Firstly,  the  average  over  a  typical  rectangle  has  Fourier  Transform: 


paJ®  =  mpw 


where,  in  the  case  of  simple  averaging,  B  is  the  scalar: 


(3.13) 


m  = 


.  UtAz  u  Ay 

4sm  — - —  sin  — — 
2  2 


u^jyAxAy 


(3.14) 


Next,  assuming  measurements: 


ffl  =  fl[3)t/0p)  4-  W  , 


(3.15) 


where  IK  is  a  white-noise  measurement  error,  the  best  estimate  p  of  pAV  is 
given,  as  in  reference  [1],  by 


where 


Hty  — 


p'  =  0  rPM<3)  . 


*ivlH{^)<t>u0(V)CT(X)BT(Z3) 
1  +<t>up)HT(12)*rfH(-V) 


(3.16) 


(3.17) 


and  and  <t>uo  are  the  spectral  densities  of  the  measurement  error  and  the 

disturbing  potential  at  the  earth's  surface. 

Combining  (3.9),  (3.11),  (3.12),  (3.15),  (3.16)  and  (3.17)  we  can  deduce  the 
spectral  density  of  the  error  A t].  After  some  simplification  (see  Appendix  A)  this 


(3.18) 


/r(Q)q-j)0L-orj)Crp)/[3) 


0yop)//rP)»^/^-a) 

1  +  tuj[z))HTia)*jjHn-3) 


E'  /rP»m-2)5(-3)q-3)0Lro(u) 


•  C  rP)B  r(3)AuKJ.m)  +  /  r(-3)[q^Hq-3)*c,p)C  rp)[fl  rp)-7]/-) 
where  V]'  denotes  2  E  *  aQd  /  denotes  the  identity  matrix. 

m,n  -  m  n 

m^O  or  VO 

The  first  term  on  the  righthand  side  of  (3.25)  is  the  spectrum  of  the  error  in 
estimating  rj  due  to  measurement  errors  without  any  mascon  approximation. 
The  terms  in  parentheses  {  }  in  the  second  term  on  the  righthand  side  of  (3.25) 
constitute  the  error  in  r]  due  to  the  mascon  approximation,  assuming  perfect 
measurements.  The  factor  in  front  of  the  parentheses  is  the  ratio 
(signal  )/(signal  +  noise  ),  and  attenuates  at  high  frequency. 


4.  NUMERICAL  RESULTS 

The  quantity  tj  of  interest  is  taken  to  be  the  vertical  gravity  anomaly  g{:) 
at  satellite  altitude  z.  Thus,  as  in  reference  [1],  the  scalar  C  rM/M  is  just 
’  where  w  =  |w|  =  \J ^ i +^y” •  The  A5(j)  and  £’U(y)  spacings  are 
assumed  equal:  =  Ay  =  A.  The  estimation  of  average  densities  is  assumed 


to  be  based  on  vertical  gradiometer  measurements  in  a  polar  satellite  at  the 
altitude  r  =  180  km.  Thus,  as  in  reference  [1],  H[3)  is  the  scalar 
Independent  gradiometer  measurements  are  assumed  to  be  available  every  8  sec 
for  6  months  with  an  accuracy  of  either  0.1  Eotvos  or  0.01  Eotvos,  so  that 
<f>lv  =  2x2/?2<y£/;V,  with  N  21  1.92  X  10®,  /?( Arm)  =  earth  radius  and 

(rw  =  10  ~10  sec'2  or  10  ~11  sec'2.  The  spectrum  of  the  surface  disturbing 
potential  is  obtained,  as  in  reference  [l],  from  the  top  3  layers  of  a  model  due  to 
Heller  containing  5  layers  of  fictitious  white-noise  potential.  The  constructed 
mascon  model  is  assumed  to  be  superimposed  on  a  low-order  potential  coefficient 
model  up  to  degree  1$  =  10.  The  mean  squared  A^  is  thus  obtained,  as  in 
reference  [l],  by  integrating  (l/4?r2)  times  <$Ar)(2),  in  (3.25),  over  the  region 

-  >  w*. 

Table  1  corresponds  to  the  integral  of  the  first  term  in  (3.25),  while  Tables 
2a  and  2b  correspond  to  the  integral  of  the  second  term.  The  effect  of  crw  in 
Tables  2  was  not  noticeable  since  the  factor  (signal)/( signal  +  noise)  was 
essentially  unity  over  the  dominant  frequency  range. 

Comparison  of  Tables  1  and  2a  suggest  that  the  ratio  A /:  need  not  be 
chosen  less  than  0.5  if  <7^  is  0.1  Eotvos,  or  not  much  less  than  0.2  if  crw  is  0  01 
Eotvos.  However,  in  practice,  a  potential  coefficient  model  is  truncated  at  some 
maximum  degree  /j  and  the  consequent  truncation  error  will  often  dominate  the 

estimation  error  in  Table  1.  The  from  truncation  is  the  integral  of 

(l/4x2)  times  (v  2e  _2w:c>fo(~'))  over  the  region  -j  >  /j / /?.  The  corresponding 


Table  1  (m  gal)  From  Estimation  Error 
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Table  2a  (m  gal)  From  Mascon  Model  Error 
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Table  2b  (m  gal)  From  Mascon  Model  Error 
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Table  3  (m  gal)  From  Truncation 
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<7^j)  is  shown  in  Table  3. 

Comparison  of  Tables  2b  and  3  suggests,  for  example,  that  a  mascon  model 
with  A  =  165  km  will  yield  about  the  same  accuracy  (.78  mgal)  at  ?  =  180  km 
as  a  potential  coefficient  model  truncated  at  degree  lx  =  90.  Since  the  number  of 
potential  coefficients  ~  =  8,100,  while  the  number  of  mascons  in  2  layers 

~  2  X  4~R  2/A  2  ~  37,470,  the  mascon  model  would  appear  to  be  less  efficient 
by  a  factor  4.  However  it  is  evident  that  in  computing  g(z)  at  a  particular 
position  in  orbit  the  mascons  located  at  a  distance  greater  than  r  from  the 
subsatellite  point  can  be  omitted  if  rvxjr  is  sufficiently  large,  without  appreciable 
degradation  in  accuracy.  To  give  a  quantitative  answer  to  the  approximate 
choice  of  it  would  be  necessary  to  take  into  account  the  cross-correlation  of 
the  mascons.  This  would  lead  to  a  quadruple  integral  for  the  spectrum  of  the 
additional  error  and  a  sextuple  integral  for  the  mean-squared  additional  error 
itself! 

A  crude  criterion  for  the  choice  of  r  is  obtainable  if  we  pretend  that  the 
mascons  are  statistically  independent.  Since  the  vertical  gravity  due  to  a  mascon 
ml}  at  relative  horizontal  position  7i}  contains  the  factor  z/{  z  2+r’;) ' ,  the 
additional  error  due  to  omission  of  mascons  with  r(;  >  is  given  by: 


[ilii 
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°\g(z)IVg(z)  =  J  rz-iz'-r^-'drl  J  rz2(z2+r2)-zdr 

0 


MAX 


=  U+^.,-A2)-2 


where  ajj.)  is  the  a  priori  variance  of  0(2).  Since  <rg^  ~  9  mgal  at 
z  =  ISO  fcm,  we  can  obtain  <7^.)  ~  0.55  mgal  by  choosing  r^  =  3.92.  With 
the  spacing  A  =  150  km,  we  need  only  about  138  local  mascons  to  achieve  a 
total  error  of  about  0.78  m  gal,  equivalent  to  that  obtained  from  the  potential 
coefficient  model  with  8,100  terms. 


5.  MASCONS  WITH  IMPROVED  AVERAGING 

Instead  of  the  scalar  averaging  operator  Bp)  of  (3.21)  corresponding  to 

uniform  averaging  over  the  typical  rectangle,  it  is  possible  to  choose  Bp)  so  as 

to  minimize  the  second  term  in  (3.25).  This  can  be  carried  out  either  by 

restricting  Bp)  in  (3.20)  to  be  a  scalar,  or  by  allowing  B(Zj)  to  be  a  2X2 

matrix.  In  either  case,  the  averaging  suffers  from  the  disadvantage  of  being 

dependent  on  /p);  i.e.  in  either  case  the  computed  mascons  would  depend  on 

satellite  altitude  2.  On  the  other  hand,  the  improvements  over  simple  averaging 

are  surprising.  Table  4  shows  the  mascon  model  error  with  optimal  scalar 

averaging  in  (3.20),  and  Table  5  shows  the  results  for  optimal  matrix  averaging. 

The  improvement  in  Table  5  over  Table  4  is  presumably  due  to  taking  better 
advancage  of  Jeffreys'  dynamic  compensation. 
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APPENDIX  A 


Firstly,  the  error  ±q(7)  like  q  *( r)  is  no  longer  a  stationary  two-dimensional 
process.  Its  autocorrelation  function  d>Af)(r^)  4  E[±q{ r),±q(7+?)\.  depends  on 
the  position  7  relative  to  the  local  rectangle,  as  does  also  its  Fourier  Transform 
the  error  spectral  density.  But  we  shall  subsequently  average  over  7. 
Now  o ^  is  expressible  as 

~~  _  "F  ^ f)t)  ’  (A.l) 

where 

00  CO 

<V,  =  TT  /  /  m  J/V  .s'  l?W  ,!-</  m*-*  ,s-!/  |  it  U  ) 

^  -CO  -00 

(A. 2) 

00 

(/  /  PT(x+Z-2?‘  ,3 ,y"  Jdl''  <Y' 

-00 

and 

E[7(x-xf  ,rf~r/  )pT{z+^-jr  ,^-q-r/'  )]  =  C-p  (£+*'  -z"  ,q+r/  f  )  , 
so  that,  using  (2.9)  and  introducing  tf  =  f+z'  -z7’,  q'  =  74-^  -3/'  : 

///}//  +  *'H  V  V  T 

47!“  -00  71  ^ 

/V  Y  Y  )dxf  dr/  dx"  dtf  df  dq> 

Averaging  over  (x,y)  removes  all  terms  from  the  double  sum  except 


n,z-1'  >  +  mfy  -  £  ) 

a  1  A  y 
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n  =  0,  m  =  0.  I'sing  the  fact  that 


«r,p>  =  vr  /  /  *“’v  '<}•,(?  v  mc  « 

-OO 


we  are  left  with: 


Similarly 


<V„P)  =  fT(  ~ 


OnnP)  =  /r(  -  2)0,rP)JP) 


(A. 3) 


(A  4) 


Next 


«,yP>  =  -V//«','-"<,“',l^(///r(y  V  in*-/  .ri  I 

47r  -00  -00 


,y-y  )rf*»  <#)(//  pr(x  +  i  -  jf‘  ,y  +  r]  -  t/' ) 


D{z  +  £  -  if*  ,y  +  q  -  f  ,r/>  )difl  dr}1  )| d^dx] 
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Averaging  over  ( x,y )  removes  all  terms  from  the  quadruple  sum  except 


n!  =  -  n,  rri  =  -  m.  This  leaves: 


<WP>  =  E£/rl-  3WrrP  + 


(A.6) 


Furthermore,  from  (3.3),  (3.15)  and  (3.16): 


KP)  =  Vr(  -3)//(  -3)^'(3)Cr(D) 


*,rP)  =  q  -  3)^t;a(3)//r(3)quJ) 


<t>77&)  =  Vf(  '^)<t>upWP)^  •  (A. 8) 

Substitution  into  (A.l)  yields 

<W3)  =  E  /r(-3)0r(-3+=J#J(^-D+^m)^t,fo(D-DBm)//r(^Bm)+4>MJ 


•v(:kzb  jyp)  -  fT(-z)c\-z)+v  (2)//rp)0p)/p) 


-/r(-3)yr(-3)M-3)^rop)Cr(3)/(3)+/r(-3)q-^)dt;a(-)qr)y(-) 


Since  will  be  obtained  as  — -  J /  <t>^r](Zi)djJI<Ljr  it  is  permissible  to  replace 

4? t 

~  +  Z?nm  by  Zj  and  3  by  3  -  ZJnm  in  each  term  of  the  summation  in  the  first 
term  in  oArJ.  Finally,  substituting  for  v(3)  from  (3.17),  we  obtain  (3.18). 
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ABSTRACT 


Efficient  Gravity  Gradient  Data  Gathering 
M.  Bilello,  J.  V.  Breakwell,  D.  B.  DeBra 


We  are  interested  in  how  one  can  separate  the  variations  in  a  gravity  field  from  the 
measurement  noise  when  making  a  survey.  Given  a  survey  pattern  in  which  the  path  of  the 
instrument  crosses  itself  (as  it  does  in  a  series  of  orthogonal  tracks),  there  are  a  discrete 
number  of  instants  at  which  the  measurements  should  be  identical.  We  have  examined  a 
number  of  different  sequences  in  generating  the  survey  pattern  to  vary  the  times  at  which 
these  iden*:cal  conditions  occur.  The  conjecture  was  that  an  appropriate  choice  of  pattern 
could  take  advantage  of  the  time  characteristics  of  the  measurement  noise  in  permitting  a 
separation  of  noise  from  gravity  data.  We  show  the  results  as  a  function  of  the  correlation 
time  of  the  measurement  noise  for  a  simple  model  of  the  gravity  field.  For  noise  varying 
from  uncorrelated  to  a  correlation  time  comparable  to  the  survey  time,  the  variation  is 
approximately  10%.  Large  differences  in  accuracy  of  reconstruction  do  not  appear  likely 
since  our  results  give  variation  between  paths  of  approximately  2%  for  two  very  dissimilar 
paths  through  the  same  grid.  Thus  the  conjecture  has  not  been  borne  out. 


Efficient  Gravity  Gradient  Data  Gatherini 


Introduction 

The  modern  interest  in  measuring  gravity  gradients  began  in  the  late  1950s  motivated 
by  determining  the  vertical  in  a  satellite.  Early  papers  considering  the  analytical  aspect 
of  gradient  determination  were  followed  in  the  next  decade  by  a  number  of  innovative 
approaches  in  how  such  an  instrument  might  be  built.  The  revolution  in  gradiometry  was 
to  make  the  measurements  in  a  moving  vehicle  and/or  in  a  satellite  without  the  gravity 
needed  for  the  geophysical  pendulum  instruments.  An  instrument  developed  at  the  Bell 
Aerosystems  was  chosen  for  field  application  for  improvements  in  navigation.  This  instru¬ 
ment  has  been  very  successful  in  its  early  field  tests  and  is  in  production  for  deployment. 
As  a  result  of  this  success  for  the  navigation  mission,  the  Defense  Mapping  Agency  (DMA) 
through  the  Air  Force  Geophysics  l  aboratory  (AFGL)  began  the  modification  of  this  in¬ 
strument  for  gravity  gradient  measurements  for  gravity  survey  work.  Many  people  have 
subsequently  contributed  to  the  development  of  a  survey  plan  and  techniques  for  utilisa¬ 
tion  of  such  an  instrument.  In  this  paper  we  explore  the  possibility  that  given  an  area 
to  be  surveyed  and  a  track  spacing  that  has  been  determined  by  the  necessary  resolution 
of  gravity  data,  there  might  be  improvements  in  accuracy  depending  upon  the  form  of  a 
grid  pattern  used  in  overflying  the  area.  The  conjecture  is  based  on  the  fact  that  instru¬ 
ment  noise,  whether  described  in  the  time  domain  or  spectrally,  may  be  different  than  the 
equivalent  noise  associated  with  gravity  fields  for  a  given  velocity  of  the  vehicle  during  the 
survey.  When  a  survey  is  performed  with  a  grid  in  which  trades  cross  each  other,  there 
are  a  discrete  number  of  crossings  at  which  the  measurements  should  be  the  same  in  both 
directions.  Different  patterns  provide  a  different  distribution  in  time  of  when  these  points 
of  identical  measurement  occur.  It  is  this  variation  in  the  distribution  and  time  which 
could  make  a  difference  in  being  able  to  separate  signal  from  noise. 

Models 

As  indicated  in  the  introduction,  the  spectral  characteristics  of  the  gravity  field  and 
of  the  instrument  will  have  an  influence  on  the  separability  of  the  gravity  information 
from  the  instrument  noise.  With  the  amount  of  experimental  data  that  exists  from  the 
laboratory  and  early  field  trials,  it  would  be  possible  to  give  a  good  model  of  the  expected 
noise  from  a  gravity  gradiometer.  However,  to  investigate  the  potential  for  improvement 
one  can  start  with  a  much  simpler  model  of  the  instrument  noise  and  vary  its  parameters 
to  see  whether  or  not  significant  improvements  are  possible.  We  have  chosen  the  latter  ap¬ 
proach  to  investigate  the  feasibility  of  improvement  with  the  expectation  that  if  significant 
improvements  appear  possible  we  would  then  improve  the  model  using  available  empirical 
data. 

Spectral  Characteristics  of  the  Field 

We  have  used  a  model  of  the  gravity  gradient  field  that  allows  us  to  determine  the 
spatial  correlations  of  the  gravity  gradient.  (J.V.  Breakwell  [1]). 


Using  an  approximation  of  flat  earth,  we  can  write: 

U(S,h)  =  e~flwU(J,o)  where  U(u,o)  is  the  Fourier  transform  of  U(x,y,o),  potential  on 
the  reference  surface  of  the  earth  &ndU(u,  h),  is  the  Fourier  transform  of  U(x,  y,  /»),  gravity 
potential  at  altitude  h. 

Then  the  gravity  gradients  components  are  given  by: 


’U„(2,hY 

U„(CS,h) 

U„(£3,h) 

U,ACS,h) 

Uju,h). 


-MW 

-jV. 


U(2,o) 


where 


u  =  (w„u /g) 

UJ  =  +  wj 

From  Heller’s  model  referenced  in  [1],  we  get  the  spectral  density  of  U(x,y,o)  with 
correlation  distance  £>,: 

<t>v.{u)  =  4>u.{u)  = 

ial 

Equation  (1)  can  be  viewed  as  a  representation  of  a  linear  system  with  U(v,  o)  as  input 
and 


H(j  u>)  =  e 


—  us* 


-Juixuj 


as  the  transfer  function. 


Then  we  can  compute  the  spectral  densities  of  the  gravity  gradient  components  at 
altitude  h: 


h) 
h) 
h) 
h ) 

tu.AZ,  h) 


By  taking  the  inverse  Fourier  transform,  we  can  determine  the  auto-correlation  func¬ 
tions  for  the  gradients,  say  Su(z,  y,  h). 

Example:  Say  we  want  to  compute  Sytt(z,  y,  A)  we  have 

h)  =  =  ^^,w4e“5“('k+D,) 

••l 

then 

y,  h)  =  £  <t>i  I  /+°°  c*"‘**»'u*e-%*h+D>'>dUadut 

im  1  J 

or 

S^,.  (r,  9,  k)  =  V  f7’  r°° 

■'o 

that  is 

Si/.,(r,  A)  =  £  <^2ir  /  w&e3w(*+D,)J0(rw)</u/ 

.tt  •'o 

Where  Jo^to)  is  a  Bessel  function  of  the  first  kind  in  rw,  in  the  special  case  of  a  flight 
path  over  a  point  grid,  we  need  to  compute  Sp  —  E[apa J],  where  ap  is  the  sequence  of 

*pi 

signals  we  want  to  estimate  ap  = 

.*PN. 

Let’s  suppose  that  we  are  measuring  the  component  UtI  of  the  gradient,  then: 

E[ap,ap,]  —  Su,,(riji  fi) 

where 

«i  =11  W  II, 

is  the  distance  between  points  P,  and  Py 

Gravity  Survey 

To  perform  a  gravity  survey,  the  craft  which  carries  the  instruments  follows  a  particular 
path.  In  the  simple  case  of  a  square  survey  area,  a  possible  strategy  is  to  fly  parallel  tracks 
as  shown  in  Figure  1. 

However,  in  order  to  remove  drifts  and  red  noise  from  the  measurements,  a  better 
way  is  to  make  cross  checks,  taking  two  measurements  at  two  different  times  at  the  same 
point.  The  grid  of  Figure  2.1  is  an  example  of  this  type  of  flight.  Also  shown  is  the  time 
of  second  crossing,  Figure  2.2.,  and  the  time  between  the  two  crossings  versus  the  point 
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of  interest,  Figure  2.3.  One  can  see  that  for  the  path  of  Figure  2.2,  the  crossings  occur 
essentially  during  the  second  half  of  the  total  survey  time  T  and  that  when  they  begin  to 
occur,  it  is  in  such  a  way  that  they  are  close  to  each  other  in  space. 

In  order  to  get  a  better  time  and  space  distribution  of  second  crossings,  a  path  such 
as  the  one  shown  in  Figure  3.1  might  be  of  interest.  Here,  a  row  or  column  is  skipped  at 
each  pass,  and  the  effect  can  be  seen  in  Figures  3.2  and  3.3.  Basically,  second  crossings 
occur  earlier  and  two  consecutive  ones  are  more  likely  to  be  spread  in  time.  Another 
advantage  of  this  kind  of  path  is  the  possibility  to  continue  to  make  measurements  while 
turning  between  two  tracks.  If  one  row  or  column  (or  more)  is  skipped,  then  the  radius  of 
curvature  in  the  turning  is  bigger,  so  that  both  the  bank  angle  and  the  induced  acceleration 
are  smaller.  This  may  allow  the  instrument  platform  to  remain  in  tolerable  perturbations 
and  compensations  may  be  possible. 

In  view  of  the  disappointing  results  that  we  are  about  to  give,  we  did  not  pursue  the 
question  of  efficiency  due  to  variations  in  the  radius  of  turns,  nor  did  we  carry  the  study 
to  include  the  effect  of  mass  attraction  and  error  modeling  on  the  instrument. 

Criteria  for  Comparison 

Our  purpose  is  to  get  an  estimate  of  the  gravity  gradient  at  the  grid  points  with  the 
smallest  error-standard  deviation.  Since  all  points  are  a  priori  of  equal  importance,  we 
take  as  the  performance  criterion  the  arithmetic  mean  of  the  standard  deviation  obtained 
at  each  point,  that  is: 


»p~  =  jv  E*< 

*  i—i 

where  <r,  =  and  P,  is  the  variance  of  the  error  in  the  gravity  quantity  at  point 
i :  Pi  =  E\(sp,  -  $,,)*].  N  is  the  number  of  points  on  the  grid.  Thus,  we  will  be  considering 
as  the  best  path  the  one  that  minimises  the  criterion  Of„. 

Theory 

The  gradiometer  output  signals  consist  of  the  sum  of  a  signal  to  be  estimated  (gravity 
quantity)  and  the  noise  inherent  in  the  instrument. 

y  =  a  +  n  where  s  is  any  one  of  the  gravity  gradient  components  and  n  is  the 
instrument  noise.  If  we  take  M  measurements  at  M  different  times,  we  have  in  vector 
form: 

yi  ni 

y  =  3  +  n  where  y  =  ;  s  =  ;  ;  n  = 

.  yu  L  . 

where  n,  is  the  instrument  noise  at  time  etc.  ... 

If  the  pattern  is  a  square  grid  with  intersecting  points,  then  M  =  2 p3  where  p  is  the 
number  of  points  on  the  side  of  the  square  grid. 
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fig.  3.3 


The  sketches  shove  show  4  by  4  grids.  The  speed  of  the 
craft  is  uniform  and  the  turning  times  are  neglected.  In  Fig. 
2.2  and  3.2,  the  time  of  the  second  crossing  at  each  point 
(from  1  to  IS)  is  plotted,  while  Figs.  2.3  and  3.3  show  the 
time  between  the  two  crossings  for  each  point  (from  1  to 
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We  assume  for  simplicity  a  linear  estimate  from  the  observations: 

i  =  Ky  where  K  is  an  M  x  M  gain  matrix.  This  is  a  smoothing  formula  where  we 
use  all  the  collected  data  to  estimate  the  gravity  quantity  at  each  point. 

The  error  in  the  estimate  is  a  =  «  —  3  or  e  =  (J  —  K)a  -  Kn. 

Then  the  covariance  matrix  of  the  error,  say  P,  can  be  computed: 

t 

P  =  E\aaT\  =  (/  -  K)S{I  -  K)t  +  KNKt  -  ( J  -  K)E[anT]KT  -  K E\naT\{I  -  KT) 
where  5  =  E[aaT ]  and  N  4  £(nnr]. 

The  gravity  signal  s  and  the  instrument  noise  being  uncorrelated,  the  formula  reduces 


P  =  (/  -  K)S(I  -  K)t  +  KNK7 

Then  we  choose  the  gain  matrix  K  that  minimises  the  trace  of  P  (least  squares 
estimate)  that  is: 

d(tr(P))  =  tr[(— (/  -  K)S  +  KN)dKT  +  dK(-S{I  -  K)T  +  NKT\  =  0 

This  yields  -(/  -  K)S  +  KN  =  0  or  K  =  5(5  +  N)~  1  whenever  the  inverse  exists 
and  in  this  case  K  exists  and  is  unique.  We  remark  that  if  (5  +  N)  is  non  invertible  then 
E[yyT]  is  non  invertible.  Minimising  every  term  P„  leads  to  the  same  gain  matrix  K.  The 
linear  least  squares  estimate  is  then  deduced  i  =  5(5  +  N)~ly.  The  performance  of  the 
estimate  is  judged  upon  the  error  covariance  matrix  and  more  precisely  on  the  diagonal 
entries  of  this  matrix.  Substituting  for  K  in  the  expression  of  P,  we  get: 

P  =  5(5  +  N)~lN 


In  addition  to  the  fact  that  n  and  a  are  uncorrelated,  we  have  implicitly  assumed  that 
a  and  n  are  sero-mean  signals.  If  this  is  not  the  case,  (£(*)  /  0  and/or  £[n]  /  0  but  still 
£[nar]  =  £[anr]  =  0),  then  the  formulae  are  modified  in  such  a  way  that  we  replace  the 
random  variables  with  their  centered  counterparts,  namely: 

where  (y  =  a  +  n) 

'  i  =  E{a)  +  K{y  -  E(y)) 

K  =  5‘(5*  +  NT1 
,  P  =  5*(5*  +  N*)~lN* 


S'±E[aaT]  -  £lsI£[sT] 


JV'*i£[nnT)  -  £[»]£[nTl 

Then  for  a  particular  pattern  that  links  times  to  points,  we  associate  the  variance 
P(tntj)  with  the  point  which  is  flown  over  at  time 

However,  for  a  grid  with  crossed  points,  it  turns  out  that  it  is  never  necessary  to  take 
the  inverse  of  the  M  x  M  matrix  (5  +  N)  because  as  can  be  expected,  there  are  a  lot  of 
redundancies  in  the  matrix  P  computed  a a  P  =  5(5  +  N)~lN.  For  example,  if  at  times  j 
and  k  the  same  point  is  flown  over  (with  j  ^  k),  then  obviously  P(t,,t/)  =  P(tk,ti)  VI;  in 
particular,  P(th «,)  =  P(tk,  tk). 

We  detail  this  in  the  next  section  on  preliminary  numerical  results. 

Numerical  Results 

We  take  for  our  example  p  =  4  and  there  are  16  =  p2  points  on  the  grid  and  we  show 
first  how  to  reduce  the  size  of  the  matrix  to  be  inverted  (5  +  N)  (the  path  lasts  M  units 
of  time). 


s«,  n»,  n*i 

s,  =  and  ”,  =  ;  also  n,  =  n,  = 

.  3>  V  aPN  .  ,  n»V  .  .  npN 

where  the  subscripts  t  stand  for  time  and  p  for  points  (N  =  p2) 
then 


j  a,  =  Fa, 
\  n,  =  Fn, 


where  F  is  the  M  x  p2  matrix  that  maps  the  points  to  the  times,  i.e.,  F(i,j)  =  1  if 
point  j  is  flown  over  at  time  and  0  otherwise. 

F  is  full  rank  and  let  P,  be  the  pseudo-inverse  of  F[Fi  is  p2  x  M)  then  we  can  write 


I  S,  =  P5,PT  /  5,  =  FS'FT 
\  N,  =  FN,F*  \  N,  =  Fn'FT 


where 


5  =  £(sarl 
N  =  £(nnr) 


From  previous  results  we  had: 


P,  =  5|(5|  +  N,)-lN, 


which  yields 


P,  =  FS,FT[FS,F*  +  FN,FT}~lFN,FT 


P,  =  FS,Ft[F(S,  +  N,)FT\~lFN,FT 


FTif'(5,  +  iV,)fT]-1F  =  (5, +  *,)">. 


Then 


P,  =  +  N,)~lN,FT  =F  P'F*  where  P,  i  SP(S,  +  WP)-‘WP. 

Pp  is  a  p2  x  p2  matrix  the  diagonal  entries  of  which  are  repeated  in  the  diagonal  of  P,. 
Pr  gives  directly  the  covariance  of  the  gravity  gradient  at  the  points  of  interest. 

For  the  numerical  example,  we  chose  a  4  x  4  grid  with  two  different  paths  and  we 
wish  to  compare  the  performances  using  the  criterion  mentioned  earlier.  We  have  first  to 
define  the  covariance  matrices  N  and  S  and  to  construct  the  F  matrix  for  the  two  different 
paths. 

The  models  used  for  the  random  signals  n  and  a  are  exponentially  correlated.  That  is, 
the  entries  of  the  covariance  matrix  Nt  vary  as  the  exponential  of  the  time  difference  and 
the  entries  of  the  covariance  matrix  S,  vary  as  the  exponential  of  the  distance,  namely: 


=  «  '  and  5p(i,»  = 

where  r  and  6  are  correlation  time  and  correlation  distance,  respectively. 

These  models  do  not  claim  to  be  accurate  but  represent  only  a  first  try  to  get  numerical 
performance. 

Then  we  compute  Pf  =  SP(SP  +  Nf)~lN,  to  determine  the  variance  of  the  error 
associated  with  the  gravity  gradient  at  each  point  of  the  grid. 

For  the  two  paths,  we  plot  the  mean  of  the  standard  deviation  versus  r(r  =  0  corre¬ 
sponds  to  a  white  noise). 

Conclusion: 

The  spectral  models  of  instrument  noise  and  gravity  gradient  signal  we  used  in  our 
simulations  may  not  be  realistic  and  this  marks  the  limitation  of  the  results  we  got.  How¬ 
ever,  in  the  special  case  of  exponential  correlated  signals,  they  allow  us  to  answer  the 
question  of  the  best  path  (among  specified  ones)  according  to  the  criterion  we  defined.  In 
terms  of  times  of  second  crossing  and  times  between  crossings,  the  two  paths  chosen  for 
the  simulation  can  be  described  as  “very"  different.  Surprisingly  enough,  the  performances 
for  the  two  paths  are  close  to  each  other  for  the  range  of  correlation  times  we  have  run. 


However,  the  gap  is  getting  wider  in  favor  of  path  1  when  the  correlation  time  gets  larger 
but  the  performance  of  path  1  is  only  1.5  %  better  for  r  =  13  units  of  time.*  In  these 
conditions,  the  choice  of  a  “better*  path  appears  not  to  be  an  issue. 

Our  final  remark  concerns  the  nature  of  the  instrument  noise.  The  way  it  has  been 
modelled  assumed  that  it  was  stationary  (in  particular  constant  variance  at  any  time);  if 
this  is  not  the  case,  quite  different  results  may  occur;  for  example,  the  importance  of  early 
crossings  increases. 
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*  1  unit  of  time  is  the  time  required  to  fly  from  a  point  to  the  next  one. 
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versus  Cor  re  1  a-fc  i  on  T  I  rne  of  iNs-fc 
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